;FUNCTION medsigma, inarr, mask, output_mask, sigma=sigma, dark=dark, average=average

pro medsigma_darks, inarr, mask, output_mask, sigma=sigma, dark=dark, average=average, atv = atv

;+
; NAME:
;       MEDSIGMA_DARKS
; PURPOSE:
;       Compute the median at each pixel across a set of 2-d images, sigma pixels
;       with a value higher than a certain sigma and recalculate the median
; EXPLANATION:
;       Each pixel in the output array contains  the median of the
;       corresponding pixels in the input arrays.   Useful, for example to
;       combine a stack of CCD images, while removing cosmic ray hits.
;
;       This routine became partially obsolete in V5.6 with the introduction
;       of the DIMENSION keyword to the intrinsic MEDIAN() function.   However,
;       it is  still useful if a input mask is needed or if sigma sigmaion is required
;       (though it is much faster to set invalid pixels to NaN values.)
; CALLING SEQUENCE:
;       outarr = MEDSIGMA([inarr, mask, output_mask, sigma=sigma ])
; INPUTS:
;       inarr  -- A three dimensional array containing the input arrays to
;                 combine together.  Each of the input arrays must be two
;                 dimensional and must have the same dimensions.  These arrays
;                 should then be stacked together into a single 3-D array,
;                 creating INARR.  INARR can also be a list of 2-D images
;                 from which the routine will create the the 3-D array.
;                 If INARR is not provided, the user is prompted to select
;                 a list of files to use.
;
; OPTIONAL INPUT:
;       mask   -- Same structure as inarr, byte array with 1b where
;                 pixels are to be included, 0b where they are to be
;                 excluded.    For floating point images, it is much faster to
;                 set masked pixels in inarr equal to !VALUES.F_NAN (see below),
;                 rather than use the mask parameter.
;
;       sigma  -- A sigma value. Any one pixel in the set that has a value
;                 greater (or smaller) than the median+sigma, will be
;                 removed from the final median calculation. The meaning of sigma
;                 changes when a dark image is provided.
;
;       dark   -- If provided, this 2-D array will be subtracted from all the
;                 data files before the median is applied. The sigma then becomes
;                 the number of photon_noise values above the median to use
;                 to reject pixels (rejection_value = sigma * SQRT(ADU-DARK)).
;
;       average-- If this flag is set, the final image will be the AVERAGE of
;                 all the included pixels instead of the MEDIAN.
;
; OUTPUTS:
;       outarr -- The output array.  It will have dimensions equal to the
;                 first two dimensions of the input array.
;
; OPTIONAL OUPUT:
;       output_mask -- Same structure as outarr, byte array with 1b
;                      pixels are valid, 0b where all the input pixels
;                      have been masked out.
;
; RESTRICTIONS:
;        Prior to V5.6, this procedure was *SLOW* because it had to loop over
;        each pixel of the image.   See notes below about an alternative with
;        CALL_EXTERNAL.
;
; EXAMPLE:
;       Suppose one wants to combine three floating point 1024 x 1024 bias
;       frames which have been read into the IDL variables im1,im2,im3
;
;       IDL> bigim = fltarr(1024,1024,3)        ;Create big array to hold images
;       IDL> bigim(0,0,0) = im1 & bigim(0,0,1) = im2 & bigim(0,0,2) = im2
;       IDL> avgbias = medsigma(bigim)
;
;       The variable avgbias will be the desired 1024x 1024 float image.
; PROCEDURE:
;       A scalar median function over the third dimension is looped over
;       each pixel of the first two dimensions.   The /EVEN keyword is used
;       with MEDIAN (which averages the two middle values), since this avoids
;       biasing the output for an even number of images.
;
;       Any values set to NAN (not a number) are ignored when computing the
;       median.    If all values for a pixel location are NAN, then the median
;       is also returned as NAN.
;
; MODIFICATION HISTORY:
;       Written by Michael R. Greason, STX, 12 June 1990.
;       Don't use MEDIAN function for even number of images.
;          W. Landsman Sep 1996
;       Mask added.  RS Hill, HSTX, 13 Mar. 1997
;       Converted to IDL V5.0   W. Landsman   September 1997
;       Use /EVEN keyword to MEDIAN    W. Landsman  September 1997
;       Rearranged code for faster execution   W. Landsman January 1998
;       Faster execution for odd number of images   W. Landsman July 2000
;       V5.4 fix for change in SIZE() definition of undefined variable
;                W. Landsman/E. Young   May 2001
;       Use MEDIAN(/DIMEN) for V5.6 or later   W. Landsman   November 2002
;       Use keyword_set() instead of ARG_present() to test for presence of mask
;           parameter  D. Hanish/W. Landsman   June 2003
;       Modified from Astron Library medarr function to add the sigma keyword
;           S. Beland December 2004
;-
; On_error,2
;                       Check parameters.

     ; check if inarr was provided
     if size(inarr,/TNAME) eq "UNDEFINED" then begin
        ; get a list of files to process
        path=cedar_getenv('CEDAR_FITS_ROOT')
        inarr = dialog_pickfile(path=path,/multi,filter="*.fits", $
                            title = 'FITS files to be medianed')
        ;if inarr[0] eq '' then return,-1
        if inarr[0] eq '' then print, "no input array"
     endif


;****
;make filename of new median'd image
;pick file base name to be used
newnames = inarr(0)
;newnames = dialog_pickfile(filter = '*.fits', $
 ;                           title = 'FITS root file to be new name base')

;amy's modifications
;strip directory path from image newname
   file_nopath0 = reverse(strsplit(newnames, '/', /extract))

file_nopath1=strjoin(file_nopath0, ' ')
	;recreating original file_nopath
	file_nopath=file_nopath0[0]   
file_begin=strsplit(file_nopath, '.fits', /extract, /regex)
;strip off .fits suffix-old way
;   file_begin = strsplit(file_nopath, '.fits', /extract, /regex)


;amy's modifications
;add write path and _ms.fits suffix
      ms_filename = file_begin + '_ms.fits'

;*****



     ; check if inarr is a string array
     if size(inarr,/TNAME) eq "STRING" then begin
        inlist=inarr
        narr = n_elements(inarr)
        ; get the size of the array from the first image
        img = mrdfits(inarr[0],0,h1,/unsigned)
        s = size(img)
        ncol = s[1]
        nrow = s[2]
        type = s[s[0]+1]
        tempimg = fltarr(ncol,nrow,narr,/NOZERO)
        tempimg[*,*,0] = img
        img=0b
        for fn=1,narr-1 do begin
           tempimg[*,*,fn] = mrdfits(inarr[fn],0,h1,/unsigned)
        endfor
        inarr = temporary(tempimg)
     endif

     s = size(inarr)
     if s[0] NE 3 then $                    ; Input array size.
            message, "Input array must have 3 dimensions"

     if !VERSION.RELEASE GE '5.6' and (N_elements(mask) EQ 0) and n_elements(sigma) eq 0 then begin
            outarr = median(inarr,dimension=3)
if keyword_set(atv) then atv, outarr
            ;return, outarr
     endif

    ;                       Create the output array.
     ncol = s[1]
     nrow = s[2]
     narr = s[3]
     type = s[s[0] + 1]
     outarr = make_array( dimen = [ncol,nrow], /NOZERO, TYPE = type )
     help,inarr
     if N_params() GT 2 then $
            output_mask = make_array (dimen = [ncol,nrow], VALUE = 1b)
     even = (narr mod 2) EQ 0

    ;                       Combine the input arrays into the output array.

     mask_given = 0b
     if keyword_set(mask) then begin
        sm = size(mask)
        if N_elements(mask) LT 4 then $
               message,'Input mask not valid... must have 3 dimensions'
        w = where(sm[0:3] eq s[0:3], cw)
        if cw eq 4 then begin
           mask_given = 1b
        endif else begin
           message,'Mask not valid... must be same shape as input cube.'
        endelse
     endif

     ; check for a dark filename or array
     dodark=0b
     if size(dark,/TNAME) ne "UNDEFINED" then begin
        if size(dark,/TNAME) eq "STRING" then begin
           ; a filename was passed -> open it
           dark = mrdfits(dark,0,h,/unsigned)
        endif
        sd=size(dark)
        if sd[1] ne ncol or sd[2] ne nrow then begin
           print,"Dimensions of dark image is different to the data -> skip dark subtraction"
           dodark=0
        endif else begin
           for i=0,narr-1 do inarr[*,*,i] -= dark
           dodark=1b
        endelse
     endif

     if n_elements(sigma) eq 0 then sigma=0
     if n_elements(average) eq 0 then average=0
     rejected=0L
     if not mask_given then begin
         ; If the /EVEN keyword is not needed, then it is faster not to use it

         if even then begin
             for j = 0l, nrow-1 do begin
                 for i = 0l, ncol-1 do begin
                     if sigma eq 0 then begin
                         outarr[i,j] = median(inarr[i,j,*],/EVEN)
                     endif else begin
                         med = median(inarr[i,j,*],/EVEN)
                         std = stddev(inarr[i,j,*])
                         if dodark then begin
                             photon_noise = sqrt(med)
                             pos = where(abs(inarr[i,j,*]-med) lt photon_noise*sigma, count)
                         endif else begin
                             pos = where(abs(inarr[i,j,*]-med) lt std*sigma, count)
                         endelse
                         if count eq 0 then begin
                             outarr[i,j] = -999.0
                         endif else if count eq narr then begin
                             if average eq 0 then outarr[i,j]=med else outarr[i,j]=avg(inarr[i,j,*])
                         endif else begin
                             rejected+=1L
                             if average eq 0 then begin
                                outarr[i,j]=median(inarr[i,j,pos],/EVEN)
                             endif else begin
                                outarr[i,j]=avg(inarr[i,j,pos])
                             endelse
                         endelse
                     endelse
                 endfor
             endfor
          endif else begin
               for j = 0l, nrow-1 do begin
                 for i = 0l, ncol-1 do begin
                     if sigma eq 0 then begin
                         outarr[i,j] = median(inarr[i,j,*])
                     endif else begin
                         med = median(inarr[i,j,*])
                         std = stddev(inarr[i,j,*])
                         if dodark then begin
                             photon_noise = sqrt(med)
                             pos = where(abs(inarr[i,j,*]-med) lt photon_noise*sigma, count)
                         endif else begin
                             pos = where(abs(inarr[i,j,*]-med) lt std*sigma, count)
                         endelse
                         if count eq 0 then begin
                             outarr[i,j] = -999.0
                         endif else if count eq narr then begin
                             if average eq 0 then outarr[i,j]=med else outarr[i,j]=avg(inarr[i,j,*])
                         endif else begin
                             rejected+=1L
                             if average eq 0 then begin
                                outarr[i,j]=median(inarr[i,j,pos])
                             endif else begin
                                outarr[i,j]=avg(inarr[i,j,pos])
                             endelse
                         endelse
                     endelse
                  endfor
             endfor
         endelse

     endif else begin

         for j = 0l, (nrow-1) do begin
             for i = 0l, (ncol-1) do begin
                 good_pixels = 1b
                 wmask = where(mask[i,j,*],cwm)
                 if cwm gt 0 then begin
                     marr = inarr[i,j,wmask]
                 endif else begin
                     good_pixels = 0b
                     output_mask[i,j] = 0b
                 endelse

                 if good_pixels then begin
                     if sigma ne 0 then begin
                         outarr[i,j] = median(marr,/EVEN)
                     endif else begin
                         med = median(marr,/EVEN)
                         std = stddev(marr) * sigma
                         if dodark then begin
                             photon_noise = sqrt(med)
                             pos = where(abs(marr-med) lt photon_noise*sigma, count)
                         endif else begin
                             pos = where(abs(marr-med) lt std*sigma, count)
                         endelse
                         if count eq 0 then begin
                             outarr[i,j] = -999.0
                         endif else if count eq narr then begin
                             if average eq 0 then outarr[i,j]=med else outarr[i,j]=avg(marr)
                         endif else begin
                             rejected+=1L
                             if average eq 0 then begin
                                outarr[i,j]=median(marr[pos])
                             endif else begin
                                outarr[i,j]=avg(marr[pos])
                             endelse
                         endelse
                     endelse
                 endif

             endfor
         endfor
     endelse

     if size(inlist,/TNAME) ne "UNDEFINED" then inarr=inlist
     print,rejected,ncol*nrow*narr,format='("Rejected ",I,"  out of ",I," pixels")'


;****

;write dark-subtracted image and original image header to file
   print, 'WRITING TO ' + ms_filename
   mwrfits, outarr, ms_filename, header
   ;mwrfits, outarr+dark, ms_filename, header

;****

if keyword_set(atv) then atv, outarr
     ;return, outarr
print, 'end'
 end
